Catalyst.jl API
Reaction network generation and representation
Catalyst provides the @reaction_network macro for generating a complete network, stored as a ReactionSystem, which in turn is composed of Reactions. ReactionSystems can be converted to other ModelingToolkit.AbstractSystems, including a ModelingToolkit.ODESystem, ModelingToolkit.SDESystem, or ModelingToolkit.JumpSystem.
When using the @reaction_network macro, Catalyst will automatically attempt to detect what is a species and what is a parameter. Everything that appear as a substrate or product in some reaction will be treated as a species, while all remaining symbols will be considered parameters (corresponding to those symbols that only appear within rate expressions and/or as stoichiometric coefficients). I.e. in
rn = @reaction_network begin
k*X, Y --> W
endY and W will all be classified as chemical species, while k and X will be classified as parameters.
The ReactionSystem generated by the @reaction_network macro is a ModelingToolkit.AbstractSystem that symbolically represents a system of chemical reactions. In some cases it can be convenient to bypass the macro and directly generate a collection of symbolic Reactions and a corresponding ReactionSystem encapsulating them. Below we illustrate with a simple SIR example how a system can be directly constructed, and demonstrate how to then generate from the ReactionSystem and solve corresponding chemical reaction ODE models, chemical Langevin equation SDE models, and stochastic chemical kinetics jump process models.
using Catalyst, DifferentialEquations, Plots
@parameters β γ
@variables t
@species S(t) I(t) R(t)
rxs = [Reaction(β, [S,I], [I], [1,1], [2])
Reaction(γ, [I], [R])]
@named rs = ReactionSystem(rxs, t)
u₀map = [S => 999.0, I => 1.0, R => 0.0]
parammap = [β => 1/10000, γ => 0.01]
tspan = (0.0, 250.0)
# solve as ODEs
odesys = convert(ODESystem, rs)
oprob = ODEProblem(odesys, u₀map, tspan, parammap)
sol = solve(oprob, Tsit5())
p1 = plot(sol, title = "ODE")
# solve as SDEs
sdesys = convert(SDESystem, rs)
sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
sol = solve(sprob, EM(), dt=.01)
p2 = plot(sol, title = "SDE")
# solve as jump process
jumpsys = convert(JumpSystem, rs)
u₀map = [S => 999, I => 1, R => 0]
dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
jprob = JumpProblem(jumpsys, dprob, Direct())
sol = solve(jprob, SSAStepper())
p3 = plot(sol, title = "jump")
plot(p1, p2, p3; layout = (3,1))
Catalyst.@reaction_network — Macro@reaction_networkGenerates a ReactionSystem that encodes a chemical reaction network.
See The Reaction DSL documentation for details on parameters to the macro.
Examples:
# a basic SIR model, with name SIR
sir_model = @reaction_network SIR begin
c1, s + i --> 2i
c2, i --> r
end
# a basic SIR model, with random generated name
sir_model = @reaction_network begin
c1, s + i --> 2i
c2, i --> r
end
# an empty network with name empty
emptyrn = @reaction_network empty
# an empty network with random generated name
emptyrn = @reaction_networkCatalyst.make_empty_network — Functionmake_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))Construct an empty ReactionSystem. iv is the independent variable, usually time, and name is the name to give the ReactionSystem.
Catalyst.@reaction — Macro@reactionGenerates a single Reaction object.
Examples:
rx = @reaction k*v, A + B --> C + D
# is equivalent to
@parameters k v
@variables t
@species A(t) B(t) C(t) D(t)
rx == Reaction(k*v, [A,B], [C,D])Here k and v will be parameters and A, B, C and D will be variables. Interpolation of existing parameters/variables also works
@parameters k b
@variables t
@species A(t)
ex = k*A^2 + t
rx = @reaction b*$ex*$A, $A --> CNotes:
- Any symbols arising in the rate expression that aren't interpolated are treated as parameters. In the reaction part (
α*A + B --> C + D), coefficients are treated as parameters, e.g.α, and rightmost symbols as species, e.g.A,B,C,D. - Works with any single arrow types supported by
@reaction_network. - Interpolation of Julia variables into the macro works similar to the
@reaction_networkmacro. See The Reaction DSL tutorial for more details.
Catalyst.Reaction — Typestruct Reaction{S, T}One chemical reaction.
Fields
rate: The rate function (excluding mass action terms).substrates: Reaction substrates.products: Reaction products.substoich: The stoichiometric coefficients of the reactants.prodstoich: The stoichiometric coefficients of the products.netstoich: The net stoichiometric coefficients of all species changed by the reaction.only_use_rate:false(default) ifrateshould be multiplied by mass action terms to give the rate law.trueifraterepresents the full reaction rate law.
Examples
using Catalyst
@parameters k[1:20]
@variables t
@species A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]), # 0 -> A
Reaction(k[2], [B], nothing), # B -> 0
Reaction(k[3],[A],[C]), # A -> C
Reaction(k[4], [C], [A,B]), # C -> A + B
Reaction(k[5], [C], [A], [1], [2]), # C -> A + A
Reaction(k[6], [A,B], [C]), # A + B -> C
Reaction(k[7], [B], [A], [2], [1]), # 2B -> A
Reaction(k[8], [A,B], [A,C]), # A + B -> A + C
Reaction(k[9], [A,B], [C,D]), # A + B -> C + D
Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]), # A+3B+4C -> 2C + 3D
Reaction(k[13], [A,B], nothing, [3,1], nothing), # 3A+B -> 0
Reaction(k[14], nothing, [A], nothing, [2]), # 0 -> 2A
Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
Reaction(k[16], [A], [B]; only_use_rate=true), # A -> B with custom rate.
Reaction(k[17]*A*exp(B), [C], [D], [2], [1]), # 2C -> D with non constant rate.
Reaction(k[18]*B, nothing, [B], nothing, [2]), # 0 -> 2B with non constant rate.
Reaction(k[19]*t, [A], [B]), # A -> B with non constant rate.
Reaction(k[20]*t*A, [B,C], [D],[2,1],[2]) # 2A +B -> 2C with non constant rate.
]Notes:
nothingcan be used to indicate a reaction that has no reactants or no products. In this case the corresponding stoichiometry vector should also be set tonothing.- The three-argument form assumes all reactant and product stoichiometric coefficients are one.
Catalyst.ReactionSystem — Typestruct ReactionSystem{V<:Catalyst.NetworkProperties} <: AbstractTimeDependentSystemA system of chemical reactions.
Fields
eqs: The equations (reactions and algebraic/differential) defining the system.rxs: The Reactions defining the system.iv: Independent variable (usually time).sivs: Spatial independent variablesstates: All dependent (state) variables, species and non-species. Must not contain the independent variable.species: Dependent state variables representing speciesps: Parameter variables. Must not contain the independent variable.var_to_name: Maps Symbol to corresponding variable.observed: Equations for observed variables.name: The name of the systemsystems: Internal sub-systemsdefaults: The default values to use when initial conditions and/or parameters are not supplied inODEProblem.
connection_type: Type of the systemnetworkproperties:NetworkPropertiesobject that can be filled in by API functions. INTERNAL – not considered part of the public API.combinatoric_ratelaws: Sets whether to use combinatoric scalings in rate laws. true by default.continuous_events: continuous_events: AVector{SymbolicContinuousCallback}that model events. The integrator will use root finding to guarantee that it steps at each zero crossing.
discrete_events: discrete_events: AVector{SymbolicDiscreteCallback}that models events. Symbolic analog toSciMLBase.DiscreteCallbackthat executes an affect when a given condition is true at the end of an integration step.
complete: complete: if a modelsysis complete, thensys.xno longer performs namespacing.
Example
Continuing from the example in the Reaction definition:
# simple constructor that infers species and parameters
@named rs = ReactionSystem(rxs, t)
# allows specification of species and parameters
@named rs = ReactionSystem(rxs, t, [A,B,C,D], k)Keyword Arguments:
observed::Vector{Equation}, equations specifying observed variables.systems::Vector{AbstractSystems}, vector of sub-systems. Can beReactionSystems,ODESystems, orNonlinearSystems.name::Symbol, the name of the system (must be provided, or@namedmust be used).defaults::Dict, a dictionary mapping parameters to their default values and species to their default initial values.checks = true, boolean for whether to check units.networkproperties = NetworkProperties(), cache for network properties calculated via API functions.combinatoric_ratelaws = true, sets the default value ofcombinatoric_ratelawsused in calls toconvertor calling various problem types with theReactionSystem.balanced_bc_check = true, sets whether to check that BC species appearing in reactions are balanced (i.e appear as both a substrate and a product with the same stoichiometry).
Notes:
- ReactionSystems currently do rudimentary unit checking, requiring that all species have the same units, and all reactions have rate laws with units of (species units) / (time units). Unit checking can be disabled by passing the keyword argument
checks=false.
ModelingToolkit and Catalyst accessor functions
A ReactionSystem is an instance of a ModelingToolkit.AbstractTimeDependentSystem, and has a number of fields that can be accessed using the Catalyst API and the ModelingToolkit.jl Abstract System Interface. Below we overview these components.
There are three basic sets of convenience accessors that will return information either from a top-level system, the top-level system and all sub-systems that are also ReactionSystems (i.e. the full reaction-network), or the top-level system, all subs-systems, and all constraint systems (i.e. the full model). To retrieve info from just a base ReactionSystem rn, ignoring sub-systems of rn, one can use the ModelingToolkit accessors (these provide direct access to the corresponding internal fields of the ReactionSystem)
ModelingToolkit.get_states(rn)is a vector that collects all the species defined withinrn, ordered by species and then non-species variables.Catalyst.get_species(rn)is a vector of all the species variables in the system. The entries inget_species(rn)correspond to the firstlength(get_species(rn))components inget_states(rn).ModelingToolkit.get_ps(rn)is a vector that collects all the parameters defined within reactions inrn.ModelingToolkit.get_eqs(rn)is a vector that collects all theReactions andSymbolics.Equationdefined withinrn, ordering allReactions beforeEquations.Catalyst.get_rxs(rn)is a vector of all theReactions inrn, and corresponds to the firstlength(get_rxs(rn))entries inget_eqs(rn).ModelingToolkit.get_iv(rn)is the independent variable used in the system (usuallytto represent time).ModelingToolkit.get_systems(rn)is a vector of all sub-systems ofrn.ModelingToolkit.get_defaults(rn)is a dictionary of all the default values for parameters and species inrn.
The preceding accessors do not allocate, directly accessing internal fields of the ReactionSystem.
To retrieve information from the full reaction network represented by a system rn, which corresponds to information within both rn and all sub-systems, one can call:
ModelingToolkit.states(rn)returns all species and variables across the system, all sub-systems, and all constraint systems. Species are ordered before non-species variables instates(rn), with the firstnumspecies(rn)entires instates(rn)being the same asspecies(rn).species(rn)is a vector collecting all the chemical species within the system and any sub-systems that are alsoReactionSystems.ModelingToolkit.parameters(rn)returns all parameters across the system, all sub-systems, and all constraint systems.reactionparams(rn)is a vector of all the parameters within the system and any sub-systems that are alsoReactionSystems. These include all parameters that appear within someReaction.ModelingToolkit.equations(rn)returns allReactions and allSymbolics.Equationsdefined across the system, all sub-systems, and all constraint systems.Reactions are ordered ahead ofEquations with the firstnumreactions(rn)entries inequations(rn)being the same asreactions(rn).reactions(rn)is a vector of all theReactions within the system and any sub-systems that are alsoReactionSystems.
These accessors will generally allocate new arrays to store their output unless there are no subsystems. In the latter case the usually return the same vector as the corresponding get_* function.
Below we list the remainder of the Catalyst API accessor functions mentioned above.
Basic system properties
See Programmatic Construction of Symbolic Reaction Systems for examples and ModelingToolkit and Catalyst Accessor Functions for more details on the basic accessor functions.
Catalyst.species — Functionspecies(network)Given a ReactionSystem, return a vector of all species defined in the system and any subsystems that are of type ReactionSystem. To get the species and non-species variables in the system and all subsystems, including non-ReactionSystem subsystems, uses states(network).
Notes:
- If
ModelingToolkit.get_systems(network)is non-empty will allocate.
Catalyst.nonspecies — Functionnonspecies(network)Return the non-species variables within the network, i.e. those states for which isspecies == false.
Notes:
- Allocates a new array to store the non-species variables.
Catalyst.reactionparams — Functionreactionparams(network)Given a ReactionSystem, return a vector of all parameters defined within the system and any subsystems that are of type ReactionSystem. To get the parameters in the system and all subsystems, including non-ReactionSystem subsystems, use parameters(network).
Notes:
- Allocates and has to calculate these dynamically by comparison for each reaction.
Catalyst.reactions — Functionreactions(network)Given a ReactionSystem, return a vector of all Reactions in the system.
Notes:
- If
ModelingToolkit.get_systems(network)is not empty, will allocate.
Catalyst.numspecies — Functionnumspecies(network)Return the total number of species within the given ReactionSystem and subsystems that are ReactionSystems.
Catalyst.numparams — Functionnumparams(network)Return the total number of parameters within the given system and all subsystems.
Catalyst.numreactions — Functionnumreactions(network)Return the total number of reactions within the given ReactionSystem and subsystems that are ReactionSystems.
Catalyst.numreactionparams — Functionnumreactionparams(network)Return the total number of parameters within the given ReactionSystem and subsystems that are ReactionSystems.
Notes
- If there are no subsystems this will be fast.
- As this calls
reactionparams, it can be slow and will allocate if there are any subsystems.
Catalyst.speciesmap — Functionspeciesmap(network)Given a ReactionSystem, return a Dictionary mapping species that participate in Reactions to their index within species(network).
Catalyst.paramsmap — Functionparamsmap(network)Given a ReactionSystem, return a Dictionary mapping from all parameters that appear within the system to their index within parameters(network).
Catalyst.reactionparamsmap — Functionreactionparamsmap(network)Given a ReactionSystem, return a Dictionary mapping from parameters that appear within Reactions to their index within reactionparams(network).
Catalyst.isspecies — Functionisspecies(s)Tests if the given symbolic variable corresponds to a chemical species.
Catalyst.isconstant — FunctionCatalyst.isconstant(s)Tests if the given symbolic variable corresponds to a constant species.
Catalyst.isbc — FunctionCatalyst.isbc(s)Tests if the given symbolic variable corresponds to a boundary condition species.
Basic reaction properties
Catalyst.ismassaction — Functionismassaction(rx, rs; rxvars = get_variables(rx.rate),
haveivdep = nothing,
stateset = Set(states(rs)),
ivset = nothing)True if a given reaction is of mass action form, i.e. rx.rate does not depend on any chemical species that correspond to states of the system, and does not depend explicitly on the independent variable (usually time).
Arguments
rx, theReaction.rs, aReactionSystemcontaining the reaction.- Optional:
rxvars,Variables which are not inrxvarsare ignored as possible dependencies. - Optional:
haveivdep,trueif theReactionratefield explicitly depends on any independent variable (i.e. t or for spatial systems x,y,etc). If not set, will be automatically calculated. - Optional:
stateset, set of states which if the rxvars are within mean rx is non-mass action. - Optional:
ivset, aSetof the independent variables of the system. If not provided and the system is spatial, i.e.isspatial(rs) == true, it will be created with all the spatial variables and the time variable. If the rate expression contains any element ofivset, thenismassaction(rx,rs) == false. Pass a custom set to control this behavior.
Notes:
- Non-integer stoichiometry is treated as non-mass action. This includes symbolic variables/terms or floating point numbers for stoichiometric coefficients.
Catalyst.dependents — Functiondependents(rx, network)Given a Reaction and a ReactionSystem, return a vector of the non-constant species and variables the reaction rate law depends on. e.g., for
k*W, 2X + 3Y --> 5Z + W
the returned vector would be [W(t),X(t),Y(t)].
Notes:
- Allocates
- Does not check for dependents within any subsystems.
- Constant species are not considered dependents since they are internally treated as parameters.
- If the rate expression depends on a non-species state variable that will be included in the dependents, i.e. in
@parameters k @variables t V(t) @species A(t) B(t) C(t) rx = Reaction(k*V, [A, B], [C]) @named rs = ReactionSystem([rx], t) issetequal(dependents(rx, rs), [A,B,V]) == true
Catalyst.dependants — Functiondependents(rx, network)See documentation for dependents.
Catalyst.substoichmat — Functionsubstoichmat(rn; sparse=false)Returns the substrate stoichiometry matrix, $S$, with $S_{i j}$ the stoichiometric coefficient of the ith substrate within the jth reaction.
Note:
- Set sparse=true for a sparse matrix representation
- Note that constant species are not considered substrates, but just components that modify the associated rate law.
Catalyst.prodstoichmat — Functionprodstoichmat(rn; sparse=false)Returns the product stoichiometry matrix, $P$, with $P_{i j}$ the stoichiometric coefficient of the ith product within the jth reaction.
Note:
- Set sparse=true for a sparse matrix representation
- Note that constant species are not treated as products, but just components that modify the associated rate law.
Catalyst.netstoichmat — Functionnetstoichmat(rn, sparse=false)Returns the net stoichiometry matrix, $N$, with $N_{i j}$ the net stoichiometric coefficient of the ith species within the jth reaction.
Notes:
- Set sparse=true for a sparse matrix representation
- Caches the matrix internally within
rnso subsequent calls are fast. - Note that constant species are not treated as reactants, but just components that modify the associated rate law. As such they do not contribute to the net stoichiometry matrix.
Catalyst.reactionrates — Functionreactionrates(network)Given a ReactionSystem, returns a vector of the symbolic reaction rates for each reaction.
(Functions to extend or modify a network)(@id apinetworkextensionandmodification)
ReactionSystems can be programmatically extended using @add_reactions, addspecies!, addparam! and addreaction!, or using ModelingToolkit.extend and ModelingToolkit.compose.
Catalyst.@add_reactions — Macro@add_reactionsAdds the reactions declared to a preexisting ReactionSystem. Note, mutates the original network.
Notes:
- To instead generate a new network by combining two existing networks use
ModelingToolkit.extend.
Example:
rn = @reaction_network begin
@parameters G
π, 2*A --> B
end
# add this reaction into rn
@add_reactions rn begin
k*A, C --> D
endCatalyst.addspecies! — Functionaddspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)Given a ReactionSystem, add the species corresponding to the variable s to the network (if it is not already defined). Returns the integer id of the species within the system.
Notes:
disablecheckswill disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
addspecies!(network::ReactionSystem, s::Num; disablechecks=false)Given a ReactionSystem, add the species corresponding to the variable s to the network (if it is not already defined). Returns the integer id of the species within the system.
disablecheckswill disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
Catalyst.addparam! — Functionaddparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)Given a ReactionSystem, add the parameter corresponding to the variable p to the network (if it is not already defined). Returns the integer id of the parameter within the system.
disablecheckswill disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
addparam!(network::ReactionSystem, p::Num; disablechecks=false)Given a ReactionSystem, add the parameter corresponding to the variable p to the network (if it is not already defined). Returns the integer id of the parameter within the system.
disablecheckswill disable checking for whether the passed in variable is already defined, which is useful when adding many new variables to the system. Do not disable checks unless you are sure the passed in variable is a new variable, as this will potentially leave the system in an undefined state.
Catalyst.addreaction! — Functionaddreaction!(network::ReactionSystem, rx::Reaction)Add the passed in reaction to the ReactionSystem. Returns the integer id of rx in the list of Reactions within network.
Notes:
- Any new species or parameters used in
rxshould be separately added tonetworkusingaddspecies!andaddparam!.
Catalyst.setdefaults! — Functionsetdefaults!(rn, newdefs)Sets the default (initial) values of parameters and species in the ReactionSystem, rn.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
setdefaults!(sir, [:S => 999.0, :I => 1.0, :R => 1.0, :β => 1e-4, :ν => .01])
# or
@parameter β ν
@variables t
@species S(t) I(t) R(t)
setdefaults!(sir, [S => 999.0, I => 1.0, R => 0.0, β => 1e-4, ν => .01])gives initial/default values to each of S, I and β
Notes:
- Can not be used to set default values for species, variables or parameters of subsystems or constraint systems. Either set defaults for those systems directly, or
flattento collate them into one system before setting defaults. - Defaults can be specified in any iterable container of symbols to value pairs or symbolics to value pairs.
ModelingToolkit.extend — Functionextend(
sys::ModelingToolkit.AbstractSystem,
basesys::ModelingToolkit.AbstractSystem;
name,
gui_metadata
) -> ReactionSystem
extend the basesys with sys, the resulting system would inherit sys's name by default.
ModelingToolkit.extend(sys::AbstractSystem, rs::ReactionSystem; name::Symbol=nameof(sys))Extends the indicated ReactionSystem with another AbstractSystem.
Notes:
- The
AbstractSystembeing added in must be anODESystem,NonlinearSystem, orReactionSystemcurrently. - Returns a new
ReactionSystemand does not modifyrs. - By default, the new
ReactionSystemwill have the same name assys.
ModelingToolkit.compose — Functioncompose(sys, systems; name)
compose multiple systems together. The resulting system would inherit the first system's name.
ModelingToolkit.flatten — FunctionCatalyst.flatten(rs::ReactionSystem)Merges all subsystems of the given ReactionSystem up into rs.
Notes:
- Returns a new
ReactionSystemthat represents the flattened system. - All
Reactions within subsystems are namespaced and merged into the list ofReactionsofrs. The merged list is then available asreactions(rs). - All algebraic and differential equations are merged in the equations of
rs. - Currently only
ReactionSystems,NonlinearSystems andODESystems are supported as sub-systems when flattening. rs.networkpropertiesis reset upon flattening.- The default value of
combinatoric_ratelawswill be the logical or of allReactionSystems.
Base.merge! — Methodmerge!(network1::ReactionSystem, network2::ReactionSystem)Merge network2 into network1.
Notes:
- Duplicate reactions between the two networks are not filtered out.
Reactions are not deepcopied to minimize allocations, so both networks will share underlying data arrays.- Subsystems are not deepcopied between the two networks and will hence be shared.
- Returns
network1. combinatoric_ratelawsis the value ofnetwork1.
Catalyst.reorder_states! — Functionreorder_states!(rn, neworder)Given a ReactionSystem and a vector neworder, reorders the states of rn, i.e. get_states(rn), according to neworder.
Notes:
- Currently only supports
ReactionSystems without subsystems.
Network analysis and representations
Note, currently API functions for network analysis and conservation law analysis do not work with constant species (currently only generated by SBMLToolkit).
Catalyst.conservationlaws — Functionconservationlaws(netstoichmat::AbstractMatrix)::MatrixGiven the net stoichiometry matrix of a reaction system, computes a matrix of conservation laws, each represented as a row in the output.
conservationlaws(rs::ReactionSystem)Return the conservation law matrix of the given ReactionSystem, calculating it if it is not already stored within the system, or returning an alias to it.
Notes:
- The first time being called it is calculated and cached in
rn, subsequent calls should be fast.
Catalyst.conservedquantities — Functionconservedquantities(state, cons_laws)Compute conserved quantities for a system with the given conservation laws.
Catalyst.conservedequations — Functionconservedequations(rn::ReactionSystem)Calculate symbolic equations from conservation laws, writing dependent variables as functions of independent variables and the conservation law constants.
Notes:
- Caches the resulting equations in
rn, so will be fast on subsequent calls.
Examples:
rn = @reaction_network begin
k, A + B --> C
k2, C --> A + B
end
conservedequations(rn)gives
2-element Vector{Equation}:
B(t) ~ A(t) + Γ[1]
C(t) ~ Γ[2] - A(t)Catalyst.conservationlaw_constants — Functionconservationlaw_constants(rn::ReactionSystem)Calculate symbolic equations from conservation laws, writing the conservation law constants in terms of the dependent and independent variables.
Notes:
- Caches the resulting equations in
rn, so will be fast on subsequent calls.
Examples:
rn = @reaction_network begin
k, A + B --> C
k2, C --> A + B
end
conservationlaw_constants(rn)gives
2-element Vector{Equation}:
Γ[1] ~ B(t) - A(t)
Γ[2] ~ A(t) + C(t)Catalyst.ReactionComplexElement — Typestruct ReactionComplexElement{T}One reaction complex element
Fields
speciesid: The integer id of the species representing this element.speciesstoich: The stoichiometric coefficient of this species.
Catalyst.ReactionComplex — Typestruct ReactionComplex{V<:Integer} <: AbstractArray{Catalyst.ReactionComplexElement{V<:Integer}, 1}One reaction complex.
Fields
speciesids: The integer ids of all species participating in this complex.speciesstoichs: The stoichiometric coefficients of all species participating in this complex.
Catalyst.reactioncomplexmap — Functionreactioncomplexmap(rn::ReactionSystem)Find each ReactionComplex within the specified system, constructing a mapping from the complex to vectors that indicate which reactions it appears in as substrates and products.
Notes:
- Each
ReactionComplexis mapped to a vector of pairs, with each pair having the formreactionidx => ± 1, where-1indicates the complex appears as a substrate and+1as a product in the reaction with integer labelreactionidx. - Constant species are ignored as part of a complex. i.e. if species
Ais constant then the reactionA + B --> C + Dis considered to consist of the complexesBandC + D. LikewiseA --> Bwould be treated as the same as0 --> B.
Catalyst.reactioncomplexes — Functionreactioncomplexes(network::ReactionSystem; sparse=false)Calculate the reaction complexes and complex incidence matrix for the given ReactionSystem.
Notes:
- returns a pair of a vector of
ReactionComplexs and the complex incidence matrix. - An empty
ReactionComplexdenotes the null (∅) state (from reactions like ∅ -> A or A -> ∅). - Constant species are ignored in generating a reaction complex. i.e. if A is constant then A –> B consists of the complexes ∅ and B.
- The complex incidence matrix, $B$, is number of complexes by number of reactions with
\[B_{i j} = \begin{cases} -1, &\text{if the i'th complex is the substrate of the j'th reaction},\\ 1, &\text{if the i'th complex is the product of the j'th reaction},\\ 0, &\text{otherwise.} \end{cases}\]
- Set sparse=true for a sparse matrix representation of the incidence matrix
Catalyst.incidencemat — Functionincidencemat(rn::ReactionSystem; sparse=false)Calculate the incidence matrix of rn, see reactioncomplexes.
Notes:
- Is cached in
rnso that future calls, assuming the same sparsity, will also be fast.
Catalyst.complexstoichmat — Functioncomplexstoichmat(network::ReactionSystem; sparse=false)Given a ReactionSystem and vector of reaction complexes, return a matrix with positive entries of size number of species by number of complexes, where the non-zero positive entries in the kth column denote stoichiometric coefficients of the species participating in the kth reaction complex.
Notes:
- Set sparse=true for a sparse matrix representation
Catalyst.complexoutgoingmat — Functioncomplexoutgoingmat(network::ReactionSystem; sparse=false)Given a ReactionSystem and complex incidence matrix, $B$, return a matrix of size num of complexes by num of reactions that identifies substrate complexes.
Notes:
- The complex outgoing matrix, $\Delta$, is defined by
\[\Delta_{i j} = \begin{cases} = 0, &\text{if } B_{i j} = 1, \\ = B_{i j}, &\text{otherwise.} \end{cases}\]
- Set sparse=true for a sparse matrix representation
Catalyst.incidencematgraph — Functionincidencematgraph(rn::ReactionSystem)Construct a directed simple graph where nodes correspond to reaction complexes and directed edges to reactions converting between two complexes.
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
complexes,incidencemat = reactioncomplexes(sir)
incidencematgraph(sir)Catalyst.linkageclasses — Functionlinkageclasses(rn::ReactionSystem)Given the incidence graph of a reaction network, return a vector of the connected components of the graph (i.e. sub-groups of reaction complexes that are connected in the incidence graph).
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
complexes,incidencemat = reactioncomplexes(sir)
linkageclasses(sir)gives
2-element Vector{Vector{Int64}}:
[1, 2]
[3, 4]Catalyst.deficiency — Functiondeficiency(rn::ReactionSystem)Calculate the deficiency of a reaction network.
Here the deficiency, $\delta$, of a network with $n$ reaction complexes, $\ell$ linkage classes and a rank $s$ stoichiometric matrix is
\[\delta = n - \ell - s\]
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
rcs,incidencemat = reactioncomplexes(sir)
δ = deficiency(sir)Catalyst.subnetworks — Functionsubnetworks(rn::ReactionSystem)Find subnetworks corresponding to each linkage class of the reaction network.
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
complexes,incidencemat = reactioncomplexes(sir)
subnetworks(sir)Catalyst.linkagedeficiencies — Functionlinkagedeficiencies(network::ReactionSystem)Calculates the deficiency of each sub-reaction network within network.
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
rcs,incidencemat = reactioncomplexes(sir)
linkage_deficiencies = linkagedeficiencies(sir)Catalyst.isreversible — Functionisreversible(rn::ReactionSystem)Given a reaction network, returns if the network is reversible or not.
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
rcs,incidencemat = reactioncomplexes(sir)
isreversible(sir)Catalyst.isweaklyreversible — Functionisweaklyreversible(rn::ReactionSystem, subnetworks)Determine if the reaction network with the given subnetworks is weakly reversible or not.
Notes:
- Requires the
incidencematto already be cached inrnby a previous call toreactioncomplexes.
For example,
sir = @reaction_network SIR begin
β, S + I --> 2I
ν, I --> R
end
rcs,incidencemat = reactioncomplexes(sir)
subnets = subnetworks(rn)
isweaklyreversible(rn, subnets)Catalyst.reset_networkproperties! — Functionreset_networkproperties!(rn::ReactionSystem)Clears the cache of various properties (like the netstoichiometry matrix). Use if such properties need to be recalculated for some reason.
Network comparison
Base.:== — Method==(rx1::Reaction, rx2::Reaction)Tests whether two Reactions are identical.
Notes:
- Ignores the order in which stoichiometry components are listed.
- Does not currently simplify rates, so a rate of
A^2+2*A+1would be considered different than(A+1)^2.
Catalyst.isequivalent — Functionisequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true)Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystems. Ignores the names of the systems in testing equality.
Notes:
- Does not currently simplify rates, so a rate of
A^2+2*A+1would be considered different than(A+1)^2. - Does not include
defaultsin determining equality.
Base.:== — Method==(rn1::ReactionSystem, rn2::ReactionSystem)Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystems. Requires the systems to have the same names too.
Notes:
- Does not currently simplify rates, so a rate of
A^2+2*A+1would be considered different than(A+1)^2. - Does not include
defaultsin determining equality.
Network visualization
Latexify can be used to convert networks to LaTeX equations by
using Latexify
latexify(rn)An optional argument, form allows using latexify to display a reaction network's ODE (as generated by the reaction rate equation) or SDE (as generated by the chemical Langevin equation) form:
latexify(rn; form=:ode)latexify(rn; form=:sde)(As of writing this, an upstream bug causes the SDE form to be erroneously displayed as the ODE form)
If Graphviz is installed and commandline accessible, it can be used to create and save network diagrams using Graph and savegraph.
Catalyst.Graph — TypeGraph(rn::ReactionSystem)Converts a ReactionSystem into a Graphviz graph. Reactions correspond to small green circles, and species to blue circles.
Notes:
- Black arrows from species to reactions indicate reactants, and are labelled with their input stoichiometry.
- Black arrows from reactions to species indicate products, and are labelled with their output stoichiometry.
- Red arrows from species to reactions indicate that species is used within the rate expression. For example, in the reaction
k*A, B --> C, there would be a red arrow fromAto the reaction node. Ink*A, A+B --> C, there would be red and black arrows fromAto the reaction node. - Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
Catalyst.complexgraph — Functioncomplexgraph(rn::ReactionSystem; complexdata=reactioncomplexes(rn))Creates a Graphviz graph of the ReactionComplexs in rn. Reactions correspond to arrows and reaction complexes to blue circles.
Notes:
- Black arrows from complexes to complexes indicate reactions whose rate is a parameter or a
Number. i.e.k, A --> B. - Red dashed arrows from complexes to complexes indicate reactions whose rate depends on species. i.e.
k*C, A --> BforCa species. - Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
Catalyst.savegraph — Functionsavegraph(g::Graph, fname, fmt="png")Given a Graph generated by Graph, save the graph to the file with name fname and extension fmt.
Notes:
fmt="png"is the default output format.- Requires the Graphviz jll to be installed, or Graphviz to be installed and commandline accessible.
Rate laws
As the underlying ReactionSystem is comprised of ModelingToolkit expressions, one can directly access the generated rate laws, and using ModelingToolkit tooling generate functions or Julia Exprs from them.
Catalyst.oderatelaw — Functionoderatelaw(rx; combinatoric_ratelaw=true)Given a Reaction, return the symbolic reaction rate law used in generated ODEs for the reaction. Note, for a reaction defined by
k*X*Y, X+Z --> 2X + Y
the expression that is returned will be k*X(t)^2*Y(t)*Z(t). For a reaction of the form
k, 2X+3Y --> Z
the expression that is returned will be k * (X(t)^2/2) * (Y(t)^3/6).
Notes:
- Allocates
combinatoric_ratelaw=trueuses factorial scaling factors in calculating the rate law, i.e. for2S -> 0at ratekthe ratelaw would bek*S^2/2!. Ifcombinatoric_ratelaw=falsethen the ratelaw isk*S^2, i.e. the scaling factor is ignored.
Catalyst.jumpratelaw — Functionjumpratelaw(rx; combinatoric_ratelaw=true)Given a Reaction, return the symbolic reaction rate law used in generated stochastic chemical kinetics model SSAs for the reaction. Note, for a reaction defined by
k*X*Y, X+Z --> 2X + Y
the expression that is returned will be k*X^2*Y*Z. For a reaction of the form
k, 2X+3Y --> Z
the expression that is returned will be k * binomial(X,2) * binomial(Y,3).
Notes:
- Allocates
combinatoric_ratelaw=trueuses binomials in calculating the rate law, i.e. for2S -> 0at ratekthe ratelaw would bek*S*(S-1)/2. Ifcombinatoric_ratelaw=falsethen the ratelaw isk*S*(S-1), i.e. the rate law is not normalized by the scaling factor.
Catalyst.mm — Functionmm(X,v,K) = v*X / (X + K)A Michaelis-Menten rate function.
Catalyst.mmr — Functionmmr(X,v,K) = v*K / (X + K)A repressive Michaelis-Menten rate function.
Catalyst.hill — Functionhill(X,v,K,n) = v*(X^n) / (X^n + K^n)A Hill rate function.
Catalyst.hillr — Functionhillr(X,v,K,n) = v*(K^n) / (X^n + K^n)A repressive Hill rate function.
Catalyst.hillar — Functionhillar(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)An activation/repressing Hill rate function.
Transformations
Base.convert — FunctionBase.convert(::Type{<:ODESystem},rs::ReactionSystem)Convert a ReactionSystem to an ModelingToolkit.ODESystem.
Keyword args and default values:
combinatoric_ratelaws=trueuses factorial scaling factors in calculating the rate law, i.e. for2S -> 0at ratekthe ratelaw would bek*S^2/2!. Setcombinatoric_ratelaws=falsefor a ratelaw ofk*S^2, i.e. the scaling factor is ignored. Defaults to the value given when theReactionSystemwas constructed (which itself defaults to true).remove_conserved=false, if set totruewill calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)Convert a ReactionSystem to an ModelingToolkit.NonlinearSystem.
Keyword args and default values:
combinatoric_ratelaws=trueuses factorial scaling factors in calculating the rate law, i.e. for2S -> 0at ratekthe ratelaw would bek*S^2/2!. Setcombinatoric_ratelaws=falsefor a ratelaw ofk*S^2, i.e. the scaling factor is ignored. Defaults to the value given when theReactionSystemwas constructed (which itself defaults to true).remove_conserved=false, if set totruewill calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
Base.convert(::Type{<:SDESystem},rs::ReactionSystem)Convert a ReactionSystem to an ModelingToolkit.SDESystem.
Notes:
combinatoric_ratelaws=trueuses factorial scaling factors in calculating the rate law, i.e. for2S -> 0at ratekthe ratelaw would bek*S^2/2!. Setcombinatoric_ratelaws=falsefor a ratelaw ofk*S^2, i.e. the scaling factor is ignored. Defaults to the value given when theReactionSystemwas constructed (which itself defaults to true).noise_scaling=nothing::Union{Vector{Num},Num,Nothing}allows for linear scaling of the noise in the chemical Langevin equations. Ifnothingis given, the default value as in Gillespie 2000 is used. Alternatively, aNumcan be given, this is added as a parameter to the system (at the end of the parameter array). All noise terms are linearly scaled with this value. The parameter may be one already declared in theReactionSystem. Finally, aVector{Num}can be provided (the length must be equal to the number of reactions). Here the noise for each reaction is scaled by the corresponding parameter in the input vector. This input may contain repeat parameters.remove_conserved=false, if set totruewill calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.- Does not currently support
ReactionSystems that include coupled algebraic or differential equations.
Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)Convert a ReactionSystem to an ModelingToolkit.JumpSystem.
Notes:
combinatoric_ratelaws=trueuses binomials in calculating the rate law, i.e. for2S -> 0at ratekthe ratelaw would bek*S*(S-1)/2. Ifcombinatoric_ratelaws=falsethen the ratelaw isk*S*(S-1), i.e. the rate law is not normalized by the scaling factor. Defaults to the value given when theReactionSystemwas constructed (which itself defaults to true).- Does not currently support
ReactionSystems that include coupled algebraic or differential equations. - Does not currently support continuous events as these are not supported by
ModelingToolkit.JumpSystems.
ModelingToolkit.structural_simplify — Functionstructural_simplify(sys)
structural_simplify(sys, io; simplify, kwargs...)
Structurally simplify algebraic equations in a system and compute the topological sort of the observed equations. When simplify=true, the simplify function will be applied during the tearing process. It also takes kwargs allow_symbolic=false and allow_parameter=true which limits the coefficient types during tearing.
The optional argument io may take a tuple (inputs, outputs). This will convert all inputs to parameters and allow them to be unconnected, i.e., simplification will allow models where n_states = n_equations - n_inputs.
Unit validation
ModelingToolkit.validate — Methodvalidate(rx::Reaction; info::String = "")Check that all substrates and products within the given Reaction have the same units, and that the units of the reaction's rate expression are internally consistent (i.e. if the rate involves sums, each term in the sum has the same units).
ModelingToolkit.validate — Functionvalidate(rs::ReactionSystem, info::String="")Check that all species in the ReactionSystem have the same units, and that the rate laws of all reactions reduce to units of (species units) / (time units).
Notes:
- Does not check subsystems, constraint equations, or non-species variables.
Utility functions
Catalyst.symmap_to_varmap — Functionsymmap_to_varmap(sys, symmap)Given a system and map of Symbols to values, generates a map from corresponding symbolic variables/parameters to the values that can be used to pass initial conditions and parameter mappings.
For example,
sir = @reaction_network sir begin
β, S + I --> 2I
ν, I --> R
end
subsys = @reaction_network subsys begin
k, A --> B
end
@named sys = compose(sir, [subsys])gives
Model sys with 3 equations
States (5):
S(t)
I(t)
R(t)
subsys₊A(t)
subsys₊B(t)
Parameters (3):
β
ν
subsys₊kto specify initial condition and parameter mappings from symbols we can use
symmap = [:S => 1.0, :I => 1.0, :R => 1.0, :subsys₊A => 1.0, :subsys₊B => 1.0]
u0map = symmap_to_varmap(sys, symmap)
pmap = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])u0map and pmap can then be used as input to various problem types.
Notes:
- Any
Symbol,sym, withinsymmapmust be a valid field ofsys. i.e.sys.symmust be defined.